
compname = getenv('computername');
if strcmp(compname,'DESKTOP-BHR0AU7')
    datadir = 'E:\ANL\Experiments\DATA\VAST\Pupillometry\';
    figDir = 'E:\ANL\Experiments\RESULTS\VAST\FIGURES\Pupillometry\';
    Rsavedir = 'E:\ANL\Experiments\RESULTS\VAST\DATAforR\WholeTrial\';
elseif strcmp(compname,'MSI')
    % use solid state external
    datadir = 'E:\ANL\PREPROCESSED_DATA\VAST\';
    figDir = 'E:\ANL\RESULTS\VAST\FIGURES\Pupillometry\';
    Rsavedir = 'E:\ANL\RESULTS\VAST\DATAforR\WholeTrial\';
else
    error('specify directories for this machine')
end

toPlot = [0; % 1) No INT, all 4 MEM -- Trial Window
          0; % 2) AS INT, all 4 MEM -- Trial Window
          0]; % 3) AT INT, all 4 MEM -- Trial Window

pupil_fs = 500;
foldernames = dir(strcat(datadir, 'Experiment_Data'));
foldernames = extractfield(foldernames, 'name');
foldernames = foldernames(3:end);
nSs = length(foldernames);
nBlk = 12;
nTrlPerBlk = 20;
[behav_a, behav_v] = deal(cell(1,nSs));
[avg_pupil_indSs, avg_pupil_indSs_res] = deal(cell(1,nSs));
allconds = {'as_none','as_as','as_at','at_none','at_as','at_at',...
            'vs_none','vs_as','vs_at','vt_none','vt_as','vt_at'};
allconds_res = allconds(~contains(allconds, '_none'));        
ENC_sample = zeros(nSs, length(allconds)); % sample of window switch in 
                                           % participant average data
ENC_sample_res = zeros(nSs, length(allconds_res)); % for residuals                                        
ENC_lengths_all = cell(size(allconds)); % ENC win lengths on each trial
% alpha_crit = 0.05;
alpha_crit = 0.01; % alpla level for significance in perm testing
       

%% Map pupil data to conditions, average within each participant

for s = 1:nSs
    if s == nSs
        fprintf('\nProcessing subject %s of %s\n', num2str(s), num2str(nSs))
    else
        fprintf('\nProcessing subject %s of %s', num2str(s), num2str(nSs))
    end
    
    %% Preallocate
    avg_pupil_indSs{s} = cell(2,length(allconds));
    avg_pupil_indSs{s}(1,:) = allconds;
    avg_pupil_indSs_res{s} = cell(2,length(allconds_res));
    avg_pupil_indSs_res{s}(1,:) = allconds_res;
    [mem_conds_a, dist_conds_a, mem_conds_v, dist_conds_v] =...
        deal(cell(nBlk * nTrlPerBlk, 1));
    
    %% Load the current behavioral datasets
    currdir = strcat(datadir,'Experiment_Data\',foldernames{s},'\');
    currfiles = dir(currdir);
    currfiles = extractfield(currfiles,'name');
    currfiles = currfiles(3:end);
    audfile = currfiles{contains(currfiles,'a')};
    visfile = currfiles{contains(currfiles,'v')};
    temp = load(strcat(currdir, audfile), 'cfg');
    behav_a{s} = temp.cfg;
    temp = load(strcat(currdir, visfile), 'cfg');
    behav_v{s} = temp.cfg;
    
    % Get condition references
    for blk = 1:nBlk
        curr_trls = (blk-1)*nTrlPerBlk + 1 : blk*nTrlPerBlk;
        for t = curr_trls
            mem_conds_a{t} = strcat(behav_a{s}.m1{blk}, behav_a{s}.d{blk});
            dist_conds_a{t} = behav_a{s}.dist{blk};
            mem_conds_v{t} = strcat(behav_v{s}.m1{blk}, behav_v{s}.d{blk});
            dist_conds_v{t} = behav_v{s}.dist{blk};
        end
    end
    
    %% Load the pupillometry data
    currdir = strcat(datadir,'Cleaned_Pupil_Data_wholeTrial\',foldernames{s},'\');
    currfiles = dir(currdir);
    currfiles = extractfield(currfiles,'name');
    currfiles = currfiles(3:end);
    audfile = currfiles{contains(currfiles,'a_cleaned')};
    visfile = currfiles{contains(currfiles,'v_cleaned')};
    
    % ----- AUDITORY -----
    temp = load(strcat(currdir, audfile));
    if sum(contains(fieldnames(temp),'cleanedData')) < 1 ||...
       sum(contains(fieldnames(temp),'metadata')) < 1
        error('Pupillometry data does not contain the expected fields')
    end
    subjpupil = temp.cleanedData.anyl.mean;
    subjmeta = temp.metadata;
    clear temp

    % Average the auditory pupillometry traces in each condition
    audinds = find(contains(allconds, 'as_') | contains(allconds, 'at_'));
    audconds = allconds(audinds);
    for c = 1:length(audconds)
        memdist_sep = strfind(audconds{c}, '_');
        currmem = audconds{c}(1:memdist_sep - 1);
        currdist = audconds{c}(memdist_sep + 1:end);
        condInds = strcmpi(mem_conds_a, currmem) &...
            strcmpi(dist_conds_a, currdist);
        currpupil = subjpupil(condInds);
        currTrigs = subjmeta.trialTrigs(condInds);
        
        % Remove empty trials (manually or automatically rejected)
        currTrigs(cellfun(@isempty, currpupil)) = [];
        currpupil(cellfun(@isempty, currpupil)) = [];
        
        % Set triggers relative to analysis win offset:
        % step 1) set triggers relative to anyl win trig [10]
        currTrigs = cellfun(@(x) [x(:,1), x(:,2) - x(x(:,1)==10,2)],...
            currTrigs, 'UniformOutput', false);
        change_win = cell2mat(cellfun(@(x) x(x(:,1)==13,2) + 30,...
            currTrigs, 'UniformOutput', false));
        % step 2) shift triggers again re ENC offset [13] + 30 (60ms stim)
        currTrigs = cellfun(@(x) [x(:,1), x(:,2) - (x(x(:,1)==13,2) + 30)],...
            currTrigs, 'UniformOutput', false);
        
        % Define the longest ENC and shortest RET windows
        ENC_lengths = cell2mat(cellfun(@(x) abs(x(x(:,1)==10,2)),...
            currTrigs, 'UniformOutput', false));
        longest_ENC = max(ENC_lengths);
        RET_lengths = zeros(size(ENC_lengths));
        for i = 1:length(RET_lengths)
            RET_lengths(i) = length(currpupil{i}) - ENC_lengths(i);
        end
        shortest_RET = min(RET_lengths);
        
        % Store the ENC window lengths (for trial count histogram)
        ENC_lengths_all{audinds(c)} = [ENC_lengths_all{audinds(c)};...
                                       ENC_lengths];
        
        % Adjust trial lengths via nan-padding / cutting; baseline correct 
        % to start of the analysis window (start ENC)
        for i = 1:length(currpupil)
            % baseline correct currpupil
            currpupil{i} = currpupil{i} - currpupil{i}(1);
            % Prepend NaNs to ensure equal ENC lengths
            NsampsBefore = longest_ENC - ENC_lengths(i);
            currpupil{i} = [nan(1,NsampsBefore), currpupil{i}];
            % Cut the ends of traces to ensure equal retention lengths
            NextraSamps = (length(currpupil{i}) - longest_ENC) - shortest_RET;
            currpupil{i} = currpupil{i}(1:end-NextraSamps);
        end
        
        % Average, now that everything is the same length
        condMean = nanmean(cell2mat(currpupil), 1);
        avg_pupil_indSs{s}{2, audinds(c)} = condMean;
        ENC_sample(s,audinds(c)) = longest_ENC;
    end
        
    % ----- VISUAL -----
    temp = load(strcat(currdir, visfile));
    if sum(contains(fieldnames(temp),'cleanedData')) < 1 ||...
       sum(contains(fieldnames(temp),'metadata')) < 1
        error('Pupillometry data does not contain the expected fields')
    end
    subjpupil = temp.cleanedData.anyl.mean;
    subjmeta = temp.metadata;
    clear temp

    % Average the visual pupillometry traces in each condition
    visinds = find(contains(allconds, 'vs_') | contains(allconds, 'vt_'));
    visconds = allconds(visinds);
    for c = 1:length(visconds)
        memdist_sep = strfind(visconds{c}, '_');
        currmem = visconds{c}(1:memdist_sep - 1);
        currdist = visconds{c}(memdist_sep + 1:end);
        condInds = strcmpi(mem_conds_v, currmem) &...
            strcmpi(dist_conds_v, currdist);
        currpupil = subjpupil(condInds);
        currTrigs = subjmeta.trialTrigs(condInds);
        
        % Remove empty trials (manually or automatically rejected)
        currTrigs(cellfun(@isempty, currpupil)) = [];
        currpupil(cellfun(@isempty, currpupil)) = [];
        
        % Set triggers relative to analysis win offset:
        % step 1) set triggers relative to anyl win trig [10]
        currTrigs = cellfun(@(x) [x(:,1), x(:,2) - x(x(:,1)==10,2)],...
            currTrigs, 'UniformOutput', false);
        % step 2) shift triggers again re ENC offset [13] + 30 (60ms stim)
        currTrigs = cellfun(@(x) [x(:,1), x(:,2) - (x(x(:,1)==13,2) + 30)],...
            currTrigs, 'UniformOutput', false);
        
        % Define the longest ENC and shortest RET windows
        ENC_lengths = cell2mat(cellfun(@(x) abs(x(x(:,1)==10,2)),...
            currTrigs, 'UniformOutput', false));
        longest_ENC = max(ENC_lengths);
        RET_lengths = zeros(size(ENC_lengths));
        for i = 1:length(RET_lengths)
            RET_lengths(i) = length(currpupil{i}) - ENC_lengths(i);
        end
        shortest_RET = min(RET_lengths);
        
        % Store the ENC window lengths (for trial count histogram)
        ENC_lengths_all{visinds(c)} = [ENC_lengths_all{visinds(c)};...
                                       ENC_lengths];
        
        % Adjust trial lengths via nan-padding / cutting; baseline correct 
        % to start of the analysis window (start ENC)
        for i = 1:length(currpupil)
            % baseline correct currpupil
            currpupil{i} = currpupil{i} - currpupil{i}(1);
            % Prepend NaNs to ensure equal ENC lengths
            NsampsBefore = longest_ENC - ENC_lengths(i);
            currpupil{i} = [nan(1,NsampsBefore), currpupil{i}];
            % Cut the ends of traces to ensure equal retention lengths
            NextraSamps = (length(currpupil{i}) - longest_ENC) - shortest_RET;
            currpupil{i} = currpupil{i}(1:end-NextraSamps);
        end
        
        % Average, now that everything is the same length
        condMean = nanmean(cell2mat(currpupil), 1);
        avg_pupil_indSs{s}{2, visinds(c)} = condMean;
        ENC_sample(s,visinds(c)) = longest_ENC;
    end
    
    % --- compute Int - No Int difference residuals --- %
    for rc = 1:length(allconds_res)
        curr_int_loc = strcmp(avg_pupil_indSs{s}(1,:), allconds_res{rc});
        noInt_name = strcat(allconds_res{rc}(1:3),'none');
        curr_noint_loc = strcmp(avg_pupil_indSs{s}(1,:), noInt_name);
        
        curr_int = avg_pupil_indSs{s}{2, curr_int_loc};
        curr_noint = avg_pupil_indSs{s}{2, curr_noint_loc};
        
        % Trim lengths to be equal
        RET_st_int = ENC_sample(s, curr_int_loc);
        RET_st_noint = ENC_sample(s, curr_noint_loc);
        % difference in RET start sample = diff ENC lengths
        if RET_st_int > RET_st_noint  % ENC longer in Int condition
            curr_int = curr_int((RET_st_int - RET_st_noint)+1:end);
        elseif RET_st_noint > RET_st_int  % ENC longer in No-Int condition
            curr_noint = curr_noint((RET_st_noint - RET_st_int)+1:end);
        end
        % now, any remaining length difference is due to a difference in
        % the retention window
        if length(curr_int) > length(curr_noint)
            curr_int = curr_int(1:length(curr_noint));
        elseif length(curr_noint) > length(curr_int)
            curr_noint = curr_noint(1:length(curr_int));
        end
        % RET start sample of the residual is equal to the earlier of the
        % two samples from the input traces
        ENC_sample_res(s,rc) = min([RET_st_int, RET_st_noint]);
        avg_pupil_indSs_res{s}{2,rc} = curr_int - curr_noint;
    end
end

save(strcat(Rsavedir,'avgPupil_wholeT_indSs.mat'), 'avg_pupil_indSs',...
    'avg_pupil_indSs_res');
save(strcat(Rsavedir,'ENC_sample_info.mat'), 'ENC_sample',...
    'ENC_sample_res', 'ENC_lengths_all');


%% Average traces across participants

pupil_indSs_temp = cell(size(avg_pupil_indSs));
for s = 1:nSs
    pupil_indSs_temp{s} = [avg_pupil_indSs{s}(1,:);...
        cell(size(avg_pupil_indSs{s}(1,:)))];
end
avg_pupil = cell(3, length(allconds));
avg_pupil(1,:) = allconds;
timebases = cell(1, length(allconds));
for c = 1:length(allconds)
    % Collect individual subjects' data in this condition, and trim the
    % data down to equate window lengths across participants. These window
    % lengths should be more or less equal across participants already, so
    % this shouldn't introduce any major timing shifts.
    short_ENC_thisC = min(ENC_sample(:,c));
    acx_ss = cell(1, nSs);
    for s = 1:nSs
        acx_ss{s} = avg_pupil_indSs{s}{2,c};
        % Remove samples on the ENC side
        N_extra = ENC_sample(s,c) - short_ENC_thisC;
        if N_extra > 0
            acx_ss{s} = acx_ss{s}(N_extra+1:end);
        end
    end
    % Now that ENC windows are equated for length, the only remaining
    % source of variability is the RET window. Trim all subject-average
    % traces from the end to be the same overall length
    short_RET_thisC = min(cellfun(@length, acx_ss));
    for s = 1:nSs
        acx_ss{s} = acx_ss{s}(1:short_RET_thisC);
        % store length-equated traces
        pupil_indSs_temp{s}{2,c} = acx_ss{s};
    end
    % Combine traces across participants in matrix form
    allCondTraces = cell2mat(acx_ss');
    % Average across participants, compute SEM
    avg_pupil{2,c} = mean(allCondTraces, 1);
    SEMtrace = zeros(1, size(allCondTraces, 2));
    for t = 1:size(allCondTraces, 2)
        SEMtrace(t) = std(allCondTraces(:,t)) / sqrt(nSs);
    end
    avg_pupil{3,c} = SEMtrace;
    % Get the timebase for the current condition
    tb = ((1:length(avg_pupil{2,c})) - 1) .* (1/pupil_fs);
    timebases{c} = tb - tb(short_ENC_thisC); % set end ENC to time zero
end


%% Compute averages of the residual traces (similar to prev block)

pupil_indSs_temp_res = cell(size(avg_pupil_indSs_res));
for s = 1:nSs
    pupil_indSs_temp_res{s} = [avg_pupil_indSs_res{s}(1,:);...
        cell(size(avg_pupil_indSs_res{s}(1,:)))];
end
avg_pupil_res = cell(3, length(allconds_res));
avg_pupil_res(1,:) = allconds_res;
timebases_res = cell(1, length(allconds_res));
for c = 1:length(allconds_res)
    % Collect individual subjects' data in this condition, and trim the
    % data down to equate window lengths across participants. These window
    % lengths should be more or less equal across participants already, so
    % this shouldn't introduce any major timing shifts.
    short_ENC_thisC = min(ENC_sample_res(:,c));
    acx_ss = cell(1, nSs);
    for s = 1:nSs
        acx_ss{s} = avg_pupil_indSs_res{s}{2,c};
        % Remove samples on the ENC side
        N_extra = ENC_sample_res(s,c) - short_ENC_thisC;
        if N_extra > 0
            acx_ss{s} = acx_ss{s}(N_extra+1:end);
        end
    end
    % Now that ENC windows are equated for length, the only remaining
    % source of variability is the RET window. Trim all subject-average
    % traces from the end to be the same overall length
    short_RET_thisC = min(cellfun(@length, acx_ss));
    for s = 1:nSs
        acx_ss{s} = acx_ss{s}(1:short_RET_thisC);
        % store length-equated traces
        pupil_indSs_temp_res{s}{2,c} = acx_ss{s};
    end
    % Combine traces across participants in matrix form
    allCondTraces = cell2mat(acx_ss');
    % Average across participants, compute SEM
    avg_pupil_res{2,c} = mean(allCondTraces, 1);
    SEMtrace = zeros(1, size(allCondTraces, 2));
    for t = 1:size(allCondTraces, 2)
        SEMtrace(t) = std(allCondTraces(:,t)) / sqrt(nSs);
    end
    avg_pupil_res{3,c} = SEMtrace;
    % Get the timebase for the current condition
    tb = ((1:length(avg_pupil_res{2,c})) - 1) .* (1/pupil_fs);
    timebases_res{c} = tb - tb(short_ENC_thisC); % set end ENC to time zero
end

save(strcat(Rsavedir,'avgPupil_wholeT.mat'), 'avg_pupil', 'timebases',...
    'avg_pupil_res', 'timebases_res');
    

%% Permutation Testing

% As a first step, this will require that all conditions have the same
% trace lengths. By this stage, traces have already been length-normalized
% across participants *within* each condition (pupil_indSs_temp).
pupil_indSs_samelen = pupil_indSs_temp;
RET_st_eachCond = min(ENC_sample, [], 1);
shortest_ENC = min(RET_st_eachCond); % this is samp 0 in permutation testing
% Pass 1: truncate the ENC windows from the beginning
for c = 1:length(allconds)
    for s = 1:nSs
        pupil_indSs_samelen{s}{2,c} =...
            pupil_indSs_samelen{s}{2,c}((RET_st_eachCond(c) - shortest_ENC)+1:end);
    end
end
% Pass 2: Truncate the RET windows from the end
length_eachCond = zeros(size(RET_st_eachCond));
for i = 1:length(allconds)
    length_eachCond(i) = length(pupil_indSs_samelen{1}{2, i});
end
shortest_RET = min(length_eachCond);
for c = 1:length(allconds)
    for s = 1:nSs
        pupil_indSs_samelen{s}{2,c} =...
            pupil_indSs_samelen{s}{2,c}(1:shortest_RET);
    end
end
clear pupil_indSs_temp

% --- Repeat this processing for the residual traces --- %
pupil_indSs_samelen_res = pupil_indSs_temp_res;
RET_st_eachCond_res = min(ENC_sample_res, [], 1);
shortest_ENC_res = min(RET_st_eachCond_res); % this is samp 0 in permutation testing
% Pass 1: truncate the ENC windows from the beginning
for c = 1:length(allconds_res)
    for s = 1:nSs
        pupil_indSs_samelen_res{s}{2,c} =...
            pupil_indSs_samelen_res{s}{2,c}((RET_st_eachCond_res(c) - shortest_ENC_res)+1:end);
    end
end
% Pass 2: Truncate the RET windows from the end
length_eachCond = zeros(size(RET_st_eachCond_res));
for i = 1:length(allconds_res)
    length_eachCond(i) = length(pupil_indSs_samelen_res{1}{2, i});
end
shortest_RET = min(length_eachCond);
for c = 1:length(allconds_res)
    for s = 1:nSs
        pupil_indSs_samelen_res{s}{2,c} =...
            pupil_indSs_samelen_res{s}{2,c}(1:shortest_RET);
    end
end
clear pupil_indSs_temp_res

% --- First Round: No Int conditions ---
% These are done on the participant average traces (not residuals)
condgroup_none = {'at_none','as_none','vt_none','vs_none'};
comparisons_none = nchoosek(1:length(condgroup_none), 2); % all combinations of 2
n_comparisons = size(comparisons_none, 1);
perm_results_none = cell(1, n_comparisons);
for i = 1:n_comparisons
    c1_idx = strcmpi(allconds, condgroup_none{comparisons_none(i, 1)});
    c1 = cellfun(@(x) x{2, c1_idx}, pupil_indSs_samelen, 'UniformOutput', false);
    c2_idx = strcmpi(allconds, condgroup_none{comparisons_none(i, 2)});
    c2 = cellfun(@(x) x{2, c2_idx}, pupil_indSs_samelen, 'UniformOutput', false);
    [perm_results_none{i}, tb_none] = pupil_permutation(c1, c2, 5000);
    tb_none = tb_none - tb_none(shortest_ENC); % set RET start to time 0
end

% --- Second Round: Int conditions ---
% These are done on the residual traces (Int - No Int)
% AT Int first...
condgroup_at = {'at_at','as_at','vt_at','vs_at'};
comparisons_at = nchoosek(1:length(condgroup_at), 2); % all combinations of 2
n_comparisons = size(comparisons_at, 1);
perm_results_at = cell(1, n_comparisons);
for i = 1:n_comparisons
    c1_idx = strcmpi(allconds_res, condgroup_at{comparisons_at(i, 1)});
    c1 = cellfun(@(x) x{2, c1_idx}, pupil_indSs_samelen_res, 'UniformOutput', false);
    c2_idx = strcmpi(allconds_res, condgroup_at{comparisons_at(i, 2)});
    c2 = cellfun(@(x) x{2, c2_idx}, pupil_indSs_samelen_res, 'UniformOutput', false);
    [perm_results_at{i}, tb_at] = pupil_permutation(c1, c2, 5000);
    tb_at = tb_at - tb_at(shortest_ENC_res); % set RET start to time 0
end

% Now AS...
condgroup_as = {'at_as','as_as','vt_as','vs_as'};
comparisons_as = nchoosek(1:length(condgroup_as), 2); % all combinations of 2
n_comparisons = size(comparisons_as, 1);
perm_results_as = cell(1, n_comparisons);
for i = 1:n_comparisons
    c1_idx = strcmpi(allconds_res, condgroup_as{comparisons_as(i, 1)});
    c1 = cellfun(@(x) x{2, c1_idx}, pupil_indSs_samelen_res, 'UniformOutput', false);
    c2_idx = strcmpi(allconds_res, condgroup_as{comparisons_as(i, 2)});
    c2 = cellfun(@(x) x{2, c2_idx}, pupil_indSs_samelen_res, 'UniformOutput', false);
    [perm_results_as{i}, tb_as] = pupil_permutation(c1, c2, 5000);
    tb_as = tb_as - tb_as(shortest_ENC_res); % set RET start to time 0
end

% Create simplified perm_results, reducing p-value to significance vectors
% 0: p > 0.05; 1: p < 0.05
[perm_results.none, perm_results.at, perm_results.as] = deal(cell(2, n_comparisons));
for cmp = 1:n_comparisons
    % None -- tests on raw (subject average) data
    perm_results.none{1,cmp} = strcat(condgroup_none{comparisons_none(cmp,1)}, '-',...
        condgroup_none{comparisons_none(cmp,2)});
    perm_results.none{2,cmp} = p_to_signif(perm_results_none{cmp}, alpha_crit, 15);
    % AT -- tests on Int-No Int residuals
    perm_results.at{1,cmp} = strcat(condgroup_at{comparisons_at(cmp,1)}, '-',...
        condgroup_at{comparisons_at(cmp,2)});
    perm_results.at{2,cmp} = p_to_signif(perm_results_at{cmp}, alpha_crit, 15);
    % AS -- tests on Int-No Int residuals
    perm_results.as{1,cmp} = strcat(condgroup_as{comparisons_as(cmp,1)}, '-',...
        condgroup_as{comparisons_as(cmp,2)});
    perm_results.as{2,cmp} = p_to_signif(perm_results_as{cmp}, alpha_crit, 15);
end
perm_TBs.none = tb_none;
perm_TBs.at = tb_at;
perm_TBs.as = tb_as;
   
save(strcat(Rsavedir,'pupilStats_wholeT.mat'), 'perm_results',...
    'perm_TBs');


% -------------------------------------------------------------------------
%% Helper functions
% -------------------------------------------------------------------------

function signif_t = p_to_signif(x, alpha, minsamps)

x = x < alpha; % convert p-values to a logical array
valshift = x(2:end)-x(1:end-1);
startSignif = find(valshift == 1) + 1;
stopSignif = find(valshift == -1);
% Case 1: whole trace is either significant or not significant
if isempty(startSignif) && isempty(stopSignif)
    if x(1) == true % all samps signif
        startSignif = 1;
        stopSignif = length(x);
    end
else % Case 2: there is some significance somewhere
    % Check the input array for boundary oddities
    % only one signif chunk, at the beginning
    if isempty(startSignif) && ~isempty(stopSignif) 
        startSignif = [1, startSignif];
    end
    % only one signif chunk, at the end
    if isempty(stopSignif) && ~isempty(startSignif)  
        stopSignif = [stopSignif, length(x)];
    end
    % first sample is true, other true chunks also present
    if startSignif(1) > stopSignif(1) 
        startSignif = [1, startSignif];
    end
    % last sample is true, other true chunks also present
    if startSignif(end) > stopSignif(end)
        stopSignif = [stopSignif, length(x)];
    end
    
    % Remove significance chunks that don't hit the minimum number of
    % consecutive samples (minsamps)
    mark_tooshort = false(size(startSignif));
    for chk = 1:length(startSignif)
        if stopSignif(chk) - startSignif(chk) + 1 < minsamps
            mark_tooshort(chk) = true;
        end
    end
    startSignif(mark_tooshort) = [];
    stopSignif(mark_tooshort) = [];
end

signif_t = [startSignif', stopSignif'];
end